Stable bright solitons in two-component Bose-Einstein condensates 
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I. INTRODUCTION 

Binary mixtures of Bose-Einstein condensates (BECs) 
are interacting quantum systems of macroscopic scale 
which exhibit rich physics not accessible in a single- 
component degenerate quantum gas. They open up in- 
triguing possibilities for a number of important physi- 
cal applications, including quantum simulation quan- 
tum interferometry and precision measurements 0, 
Experimentally, multi-component BECs are gener- 
ated as mixtures of atoms in different hyperfine states 
or by simultaneously trapping different atomic species. 
The key difference between multi-component and single- 
component BECs is the inter-component interaction. It 
is significant that the strength and even the sign of the 
various atomic interactions can be controlled by means 
of the Feshbach resonance [1, @ . 

Various types of nonlinear matter wave structures have 
been predicted for multi-component BECs 0-[Ili- How- 
ever, many basic properties of two-component BECs re- 
main to be investigated even in simple mean-field ap- 
proaches based on coupled Gross-Pitaevskii equations 
(GPEs) with a vector order parameter. Most previ- 
ous work on two-dimensional (2D) and three-dimensional 
(3D) vector solitonic structures focused on BECs with re- 
pulsive intra-componcnt interactions [T2| - [l6| . 

As is well known, bright 2D and 3D solitons, described 
by a GPE with attractive nonlinearity, are unstable and 
may collapse in a finite amount of time, if the number 
of atoms in the condensate exceeds a critical value (see 
e.g. Ref. [13] )■ Recent theoretical investigations [l^ 
predicted the existence of stable soliton-vortex pairs in 
trapped BECs with attractive intra-componcnt interac- 
tions. However, the ground state 2D soliton-soliton struc- 
tures in such BECs have never been investigated and will 
be the focus of the present work. Note, that in a differ- 
ent context soliton-soliton pairs, described by a similar 
model, but without linear potential, are found to be un- 
stable [l^ . It is reasonable to expect that a stabilization 
of fundamental vector solitons could occur in an addi- 
tional external trapping potential. 

In the present work we perform a detailed theoreti- 
cal analysis of nonlinear matter-wave structures in bi- 



nary mixtures of atoms with attractive intra-component 
interactions and attractive as well as repulsive inter- 
component interactions. General properties of the steady 
states of such systems are investigated by means of a vari- 
ational analysis and numerical simulations. The condi- 
tions for the existence and stability of matter- wave vector 
solitons are revealed. 



II. MODEL 

We consider a binary mixture of BECs at zero tem- 
perature described in mean-field approximation by two 
coupled Gross-Pitaevskii equations (GPEs) (j G {1,2}) 



ih 



dt 
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with Hj = — 2^V^-|-V, (r), and Mj the mass of an atom 
of type j loaded into the axially-symmetric harmonic ex- 
ternal trapping potential Vj{r) = MjUj\{x^ + y^)/2 + 
MjuJlz^ /2. Interactions between atoms of the same type 
(mira-component interactions) are characterized by the 
diagonal coupling coefficients gjj = inh^ajj/Mj, where 
Ujj are the s-wave scattering lengths for binary colli- 
sions between these atoms. /n<er-componcnt interac- 
tions are controlled by the off-diagonal coupling terms 
<?12 = 521 = 27r?i^ai2/M* with the reduced mass Af* — 
M1M2/ {Ml 4-M2). Note, that our simple mean field ap- 
proach provides a reasonable approximation for the evo- 
lutionary scenarios investigated here, however, a mean- 
field approach of course fails at the final stage of the 
collapse when the atomic density increases catastrophi- 
cally. 

We assume that the longitudinal trapping frequency 
ujz is much larger than the transversal trapping fre- 
quency ujj_ {ujz 3> wj.), and that the nonlinear inter- 
actions ^ Qij are weak with respect to the confine- 
ment strength of the potential in the longitudinal direc- 
tion. In this case the BEC is "disk-shaped" and we may 
assume that the longitudinal motion of condensates is 
frozen in, ^'^(r, f) = '■^ j{x,y,t)T j{z,t), where Tj{z,t) = 
(^zjVtt)^^^^ exp(— l^xt — /llj) is the ground state 
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wave function of the longitudinal harmonic trapping po- 
tential, Izj = y/Ti/ (MjUJz). After integrating out the lon- 
gitudinal coordinates, the GP equations take the effective 
2D form 



where hi 



A 



and Aj_ 



VJ-, = 0, 



(2) 



d"^ jdx^ + jdy^ is the 



2D Laplacian. Here we have introduced the dimen- 
sionless variables {x,y) — 7> {x,y)/l^_i, t t/r, and 
*j7yc7, where = y/h/{Miui±), t = 2/ujj_, 
Cj — hu!±y^Tr{lli + II2) I (2|5jj I)- Dimensionless cou- 
pling parameters are defined by 612 

621 = ^jfSiy^ 1~Pk' simplicity, we consider here only 
the symmetric case Mi = M2 and = 1722, so that 
K = 1 and 612 = 621 = f- Furthermore, we specifi- 
cally consider only cases where the diagonal part of the 
interaction matrix gij is attractive {gn < 0, 322 < 0) 
while the off-diagonal terms 312 = 521 may be attractive 
(512 < 0, tr > 0) or repulsive {gi2 > 0, cr < 0). 

The coupled differential equations ([2]) have the follow- 
ing integrals of motion: 

(i) the number of particles in each component 



(ii) the total energy 

E = Ei+E2-cT J |*in^2prf2r, 

with 



(3) 
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(iii) momentum, and (iv) angular momentum. 



III. STATIONARY SOLUTIONS 



we then obtain a two-parameter family of vector soliton 
solutions (with the chemical potentials /ii and fi2 as pa- 
rameters). 

Obviously, vector solitons with the same chemical po- 
tential = /i2 = /i and the same radial profiles ipj{r) 
are described by a single GPE with a harmonic poten- 
tial and an effective interaction strength a ~ 1 + a. It 
is known (see e.g. Ref. [l^l) that 2D solitonic solutions 
of this equation exist for < 2 if > and for fi > 2 
if o" < 0. The value /2 = 2 coincides with the eigenvalue 
of the ground state of a linear Schrodinger equation with 
harmonic oscillator potential, i.e. nonlinear terms vanish 
and the number of atoms tends to zero when /2 — >■ 2. If 
(7 < 0, the number N of atoms grows rapidly for /i > 2; 
i{ (T > 0, N saturates at a critical value Ncr/a- where 
/2 — > — 00. Here, A^cr = 11-68 is the number of atoms in 
the fundamental Townes soliton, which is the localized 
solution of a single GPE without external trap. 

Using these results for one-component solitons as a 
guide, one may expect qualitatively different properties 
for vector solitons with strong repulsion (cr < —1) or 
weak repulsion — 1 < cr < 0. Analogously, we will also 
discuss separately the properties of the soliton-soliton 
pairs with weak attraction (0 < cr < 1) and strong at- 
traction (cr > 1). 



A. Attractive inter-component interactions (cr > 0) 

The coupled differential equations © were solved nu- 
merically by a relaxation method similar to that de- 
scribed in Ref. [2l|. Furthermore, we performed a vari- 
ational analysis of the fundamental vector solitons. The 
variational results not only provide an appropriate initial 
condition for our numerical relaxation procedure but also 
they show a good agreement with our numerical results 
as is seen e.g. from Figs. [1] (b)-(d) and (f)-(h). 

If all interactions between particles are attractive, then 
both solitonic components are basically bell-shaped and 
a simple trial function of the form 



A, 



(7) 



We consider now stationary vector soliton solutions of 
the coupled GPEs ([2]). We are looking for a ground state 
of the form 



*,(r,t)=Vj(r)e-'''^* 



(5) 



where r = y^x"^ + y"^. Each solution is characterized by 
chemical potentials ^j. The real functions ipj{r) satisfy 
the coupled equations 

l^i^^j + + ^V-j- - r"^, + {^] + V', = (6) 

and the boundary conditions V'j(O) = and i\}j{co) ~ 0. 
At fixed strength of the inter-component interaction cr 



is expected to be a good approximation. Here aj is the 
effective width of component j. The trial function is 
normalized to the number of particles: — ^j- 

Substituting the trial function Q into Eq. ([4]), we obtain 
the energy 



E — E\ + E2 — 0E12 



(8) 



with 



E, = A^,- 



E12 



N1N2 



■K{ai 
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A soliton solution corresponds to the stationary point of 
the total energy at the fixed number of particles in each 



component: (dE/dai 



'Ni,N2 



0, {dE/da2) 



0. 
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FIG. 1: (Color online) (a), (e) Existence region for vector 
solitons in the (/xi, 11-2) plane for cr = 0.5 and a = 2.0. The 
insets show examples for radial profiles ipi (r) (red dashed line) 
and ip2{r) (black solid line) at points A, B, and C in the 
(^1, 112) plane: A = (-6, -4), B = (-5, -5), C = (-2, -8). 
(b)-(d), (f)-(h) The number of atoms A'^i (red dashed line 
for numerical, red dotted line for variational results) and 
(black solid line for numerical results, black dashed-dotted 
line for variational results). The number of atoms are plotted 
as functions of the difference /i_ = /ii — fj.2 at fixed total 
chemical potential /i+ = jJ-i + ^^2■ The dotted straight lines 
in (a) and (e) correspond to the ^+ of the particle number 
profiles given in (b) -(d) and (f)-(g), respectively. 



only single-component (scalar) solitons exist. For differ- 
ent fixed values of /i+ we determined the number of par- 
ticles in each component and then constructed the com- 
plete existence region, which is shown in Figs. [T] (a) and 
(e) . In these figures the solid black line separates the vec- 
tor solitons from the scalar soliton region where one soli- 
tonic component vanishes and the red dashed line indi- 
cates the boundary where the other component vanishes. 
As expected, these curves merge at the point fii = /i2 = 2 
where both components vanish. With increasing interac- 
tion strength a the region of existence for vector solitons 
gradually shrinks, and at tr = 1 the boundaries degener- 
ate to a straight line. Indeed, only scalar solitons with 
equal chemical potential exist for (7=1. For ct > 1, 
where inter-component interactions dominate over intra- 
component interactions, the region of existence for vector 
solitons grows again. With respect to the case a < 1 the 
boundary lines of the existence region swap around [com- 
pare Fig. [1] (a) and Fig. [T] (e)]. 

As is seen from the Fig. [1] the number of particles 
Nj (/i_ ) reaches its maximum value Nj ,„ax at the bound- 
ary of the existence domain where the other component 
A^3_j vanishes. For decreasing total chemical potential 



Nj max increases towards 



11.68 for numerical 
Att for the variational trial func- 



solutions and A^cr 
tions dZl). 

The remarkably different behavior of the functions 
Nj{^-) for weak or strong inter-component interactions 
is evident from a comparison of Figs[T](d) and (h). In the 
first case (0 < cr < 1) the number of atoms saturates at 
Njmax < Acr. lu thc sccoud casc (cr > 1) the number of 
atoms decays rapidly from Nj max ^ A"cr at one boundary 
to zero at the other boundary of the existence domain. 
It is easy to verify that for the case a > I the sum of 
the number of particles in each component is always less 
than the critical value A'cr. For weak inter-component 
interactions (0 < cr < 1) we may find A^i + A^2 > A^cr- 

In summary, if inter-component as well as intra- 
component interactions are attractive, then both soli- 
tonic components are localized at the bottom of the po- 
tential trap, and the number of atoms in each component 
is subcritical. These facts change dramatically for repul- 
sive inter-component interactions, which we discuss next. 



B. Repulsive inter-component interactions (cr < 0) 



Obviously, the two-parameter family of vector soliton 
solutions is symmetric with respect to the simultaneous 
interchange /^i -O- /Lt2 and ipi ^ "02, so that in the plane of 
chemical potentials (/xi, ^2) the solutions are symmetric 
with respect to the line fJ-i = fJ-2- To utilize this symmetry 
we present the number of particles for each component 
as functions of the difference = fJ-i — fJ-2 at the fixed 
total chemical potential /i+ = fii + fi2- Typical results 
are shown in Fig. [I](b)-(d) and (f)-(h). 

Note, that the region of existence for vector soliton 
pairs has boundaries in the (/ii, (12) plane, beyond which 



To gain insight into the properties of the fundamental 
vector solitons with inter-component repulsion we first 
present the results of a variational analysis. They provide 
a rather complete overview of the steady states which 
correspond to stationary points of the total energy (|4]) at 
a fixed number of atoms in each component. 

Minima of the total energy at a given number of atoms 
can be attained for two different types of the solutions. 
The first type we call "bell-shaped" since the density dis- 
tribution of both components looks Gaussian-like. To 
stabilize such a "miscible" configuration for negative cr. 
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a=-0.5 a=-2.0 




FIG. 2: (Color online) Variational results for a = —0.5 [(a)- 
(c)] and a = —2.0 [(d)-(f)]. Shown are the number of atoms 
A'^i (red dashed curve for solutions with 5 = 0, red curves with 
open circles for solution with 5 > 0) and N2 (solid black curves 
for solutions with 5 = and black curves with filled dots for 
solutions with 5 > 0) as the functions of /i_ = ^1 — ^2 at 
different values of the total chemical potential — fii + fi2- 
Green dotted line: N = Na = 47r. 



the particle number in each component must increase 
with respect to the previously discussed case (t > in 
order to compensate for the additional inter-component 
repulsion by increased intra-component attraction. So- 
lutions of the second type will be referred to as "phase- 
separated" or "immiscible" . These solutions minimize 
the "component-overlapping" part of the total energy 
[the last term in Eq. ^] by substantially changing the 
shape of the density distributions. One component is 
pushed outwards and forms a ring-like shell, while the 
other component is noticeably compressed - it has a 
higher peak density and narrower width than its non- 
interacting counterpart. The "immiscible" states are 
typical for BECs with both rep ulsive intra- and inter- 
component interactions [23 - l27| . The notable feature 
of BECs with attractive intra-component but repulsive 
inter-component interactions is the possible coexistence 
of "miscible" and "immiscible" states for the same chem- 
ical potentials. 

The simplest normalized trial function tpj (r) support- 
ing spatially separated components is given by the Ansatz 

V',M=A, (^l + <5,^je"^^, (9) 
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FIG. 3: (Color online) Variational parameter 5 as a function 
of /i_ = /ii — /i2 for different values of the total chemical 
potential = A^i + ^2. 

where aj is the effective width of compo- 
nent j of the soliton-soliton pair, and Aj = 

^Nj/{TTaj{l + 2Sj{l + Sj))). The parameter 5j > 

introduces deviations from the Gaussian-like shape of 
the soliton. li Sj > 1/2, the density distribution has a 
local minimum at the bottom of the external potential 
trap. Evidently only one of the two BEC components 
is forced outwards by repulsive inter-component inter- 
actions. Consequently, 62 may be set to zero without 
loss of generality and we need to determine only three 
variational parameters: ai, 02 and 6 = 61. Variational 
results for the number of atoms as functions of the 
chemical potential difference ^_ = fii — /i2 are shown in 
Fig.H 

Let us first discuss the properties of the "bell-shaped" 
solutions with 6 = 0. These solutions are presented in 
Fig. [2] by the red dashed curves for A^i and by the solid 
black line for N2. As it should, solutions with the same 
chemical potential and same number of atoms only exist 
for ^+ < 4 if -1 < cr < and for 4 if cr < -1. For 

— 1 < cr < 0, the existence domain of vector solitons is 
bounded in the (/xi, plane like in the case a > 0. 

However, if inter-component repulsion dominates over 
intra-component attraction (cr < —1) "bell-shaped" so- 
lutions exist for any chemical potential. It is remarkable 
that even at = the solitonic components differ es- 
sentially from each other for jj,^ < 4 [see Fig. [2] (f)], 
but for /i+ > 4 there appears an additional branch of 
solutions that share the same profile at = [see Fig. 
[5] (d),(e)]. The gap between these two branches grad- 
ually disappears when 11+ increases [compare Figs. [2] 
(e) and (d)]. For relatively weak inter-component inter- 
actions (— 1 < CT < 0) and decreasing ^+ the absolute 
value of the derivatives \dNj/d^ij\ grows rapidly in the 
vicinity of the intersection point. Furthermore, it is in- 
teresting, that in the profile for Nj{fi^) a fold appears 
for Nj{0) > 167r/(3 + a) provided that -1 < a < -1/3. 
In this case two additional crossing points are observed 
at /ii = ^2 [see Fig. [21(c)]. 

The properties of the second type of solutions ( "phase- 
separated" states with 6 > 0), predicted by the varia- 
tional method, are illustrated in Fig. [2]for Ni{fi^) by red 
lines with open circles, and for N2{p--) by black curves 
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FIG. 4: (Color online) Examples of radial profiles tpi (red 
dashed curves) and ijj2 (black solid curves) for a = —0.5 and 
a — —2.0 found numerically. 



with filled dots. Recall that we show only solutions with 
di — 6 and $2=0 but the rest of solutions can be eas- 
ily constructed using the corresponding symmetry. The 
parameter S describing the deviation of the radial profile 
from the Gaussian-like shape is given in Fig. [31 Here, 
(5(/j,_) is shown for the same values of a and /x+ as in 
Fig. [5] As is seen from Fig. [3] the variational parameter 
6 reaches zero at some value of /i_ , where the two types 
of solutions S > and (5 = merge. Indeed this is also 
observed in Fig. [2] for solutions with fii ~ fi2- 

The above described features of phase-separated 
steady-states obtained using the variational approach are 
supported by numerical simulations of Eqs. Note 
that only half of the phase-separated solutions are shown 
in Figs. [S](f)-(h) to avoid confusion. The traces in these 
figures can easily be completed using the replacements 
A^i — > N2 and /i_ — > — /i-. Unfortunately, our relax- 
ation technique is not able to reproduce both types of 
solutions predicted by variational method for repulsive 
inter-component interactions since the "bell-shaped" so- 
lutions (with 51=^2= 0) usually have higher energy 
than the phase-separated solutions at the same chemical 
potential. That is why our numerical procedure rather 
converges to the "phase-separated" solutions or to scalar 
solitons. Therefore, our analysis of the steady-states with 
(5i = (^2 = is limited to the variational method. 
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FIG. 5: (Golor online) Top: Stability region in the {fii, ^2) 
plane (green) as obtained by numerical simulations: (a) 
(7 = —0.5, (b) cr = —2.0. Bottom: Number of atoms A^i 
(red dashed curves) and (black solid curves) as a function 
of the chemical potential difference /i_ ~ — fJ^2 at fixed 
total chemical potential fi+ = /ii + ^2 obtained numerically: 
(c), (d), (e) a = -0.5, (f), (g), (h) a = -2.0. The dotted 
straight lines in (a) and (e) correspond to the fi+ of the parti- 
cle number profiles given in (b) -(d) and (f)-(g), respectively. 



Typical radial profiles found numerically are shown in 
Fig. m Generally, the numerical results are found to be 
in good agreement with the variational predictions for 
the "phase-separated" solutions. One expects, however, 
differences between numerical and variational results for 
solutions with a very dense ring-shaped component. In- 
deed, a ring component with large amplitude could push 
away the atoms at the periphery of the core compo- 
nent. Consequently, an additional outer ring should be 
observed in the density distribution of the core compo- 
nent. Moreover, if the number of atoms increases even 
more, then more and more rings should appear in both 
solitonic components. Actually, an additional outer ring 
is seen for '02(7') in Fig. [5] (f), and such a modification 



6 



of the core component explains the non-monotonic be- 
havior of the diagrams iV2(/i-) at negative /i_ in Figs. [5] 
(f)-(h)- 



IV. STABILITY ANALYSIS 

A few general conclusions concerning stability follow 
from our variational analysis. In fact, it is easy to 
find a sufficient condition for stability against radially- 
symmetric collapse: if both solitonic components satisfy 
the condition Nj < Ncr, then the stationary solutions 
are expected to be stable since they correspond to the 
minimum of the total energy. As was shown in the pre- 
vious section, for attractive inter-component interactions 
(cr > 0) the particle numbers in both solitonic compo- 
nents are always below the critical value and, therefore, 
the vector solitons should be stable in this case. 

For (T < 0, the variational analysis predicts the ex- 
istence of two types of stationary solutions: "phase- 
separated" states and colocated "bell-shaped" states 
with different stability properties. The solutions with 
Si ^ 62 ~ should be stable provided iVi < A^cr and 
-^2 < ^cr- At the same time, the phase-separated states 
are predicted to be stable against collapse even if the 
ring-component j (with Sj > 0) is over-critical provided 
that the bell-shaped component (with S^^j = 0) consists 
of an undercritical number of particles. 

It is clear that a radially-symmetric variational analy- 
sis does not provide any information about the stability 
with respect to azimuthal perturbations. This question 
may be addressed by numerical simulations as discussed 
below. In contrast to the results of the variational analy- 
sis our numerical simulations predict the phase-separated 
states to be unstable against collapse even if only one 
component has a number of particles above the critical 
value A^cr- Note that numerical time evolution of the vari- 
ational Ansatz for the steady-states with i5i = ^2 = con- 
firmed their stability only in the region where these solu- 
tions practically merge with the second solitonic branch 
[see Figs. HEl and[S]. 

The stability of the stationary solutions with respect 
to azimuthal symmetry-breaking perturbations was in- 
vestigated by a linear stability analysis. We write the 
order parameter of the two-component BEC in the form 
^j{r,t) = {ipj{r) + ej{r,t))e~'^^* with a small perturba- 
tion ej{r,t) = Uj(r)e''^*+*-^'^-hw*(r)e-*'^'*-*-^'^ and insert 
it into Eq. ([5]). We linearized the resulting equations with 
respect to e and solved the resulting eigenvalue problem 
for Lu. Solutions exhibiting an imaginary part would in- 
dicate instability. 

We find that jl = Im(a;) are identically zero for all L 
if the inter-component interaction is attractive. If it is 
repulsive the azimuthal modes L = 1 and L ~ 2 may be 
unstable. A typical example for a maximum growth rate 
7l as a function of the chemical potential /ii at fixed ^2 
is shown in Fig. [HI Note that above some threshold the 
growth rates vanish. It is important to realize that the 
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^^1 

FIG. 6: (Color online) Maximum growth rates as a function 
of fii at fixed ^2- Note that all growth rates vanish above 
some threshold value of fj,i. 



mode L = 1 has a much wider instability region than 
L = 2, and that the growth rates of the mode L = 2 
are smaller than those of the mode L = 1. Therefore, 
the L = 1 modes are most dangerously infiuencing the 
stability of repulsively interacting solitons. 

The predictions of the variational analysis and the re- 
sults of the linear stability analysis have been tested by 
numerical experiments of the time evolution of the per- 
turbed vector solitons. For the numerical solution of the 
time-dependent GPEs ^ a standard split-step Fourier 
transform method (see, e.g., Ref. [2I]) has been used. 
While vector solitons with attractive inter-component in- 
teractions are confirmed to be stable over the full ex- 
istence domain, repulsive inter-component interactions 
lead to various instabilities. Typical examples for un- 
stable dynamics are given in Fig. [71 

We found that, if at least one solitonic component con- 
tains an overcritical number of atoms, then the stronger 
component finally collapses. However, we do not see a 
collapse, if both solitonic components are under-critical. 
Obviously, the potential trap prevents the unlimited ex- 
pansion of the BEC. Furthermore, the final collapse is not 
possible, if the number of atoms in both components is 
subcritical, so as the result of symmetry-breaking insta- 
bility we see a contraction and relative motion of the two 
solitonic components. Oscillations between such a quasi- 
coUapse and an asymmetric state is shown in Fig. [7] (a) 
for times t = 3.5 and t — 4.0. Here, both solitonic com- 
ponents are under-critical, and the growth rate for the 
L ~ 1 mode is non- vanishing. 

The growth of the L = 2 mode leads to the decay 
of the shell component of the soliton into two filaments. 
At the beginning, two humps and two holes appear in 
the ring, and the inner component rapidly leaks through 
these hollows into a ring-shaped "well" . Repulsive inter- 
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FIG. 7: Time evolution of tlie density distributions I'l'i |^ (up- 
per rows) and |*I'2|^ (lower rows) of perturbed vector solitons 
for (a) (7 — —0.5, fii = —1.5, fi2 ~ 2 and (b) a = —0.8, 
m = —5, = -5. 

component interactions between the ring-shaped compo- 
nent and the leaking inner component leads to an in- 
crease of the hollows and finally the ring breaks up into 
two separate parts as seen in Fig. [7] (b) . Note that the 
dynamics of the unstable trapped two-component BEC 
is similar to the evolution of vector solitons reported 
in Rcf. [l9| within a similar model but without exter- 
nal trapping potential. The most outstanding feature 
of trapped matter-wave vector solitons is the fact that 
they can be completely stabilized even for strong repul- 
sive intcr-component interactions. We show the stabil- 
ity region in green color in the {fii, (12) plane for weak 
(a ~ —0.5) and strong (a = —2.0) intcr-component re- 



pulsive interactions in Figs. [S] (a) and (e). 



V. SUMMARY AND CONCLUSIONS 

Fundamental 2D soliton-soliton pairs arc investigated 
in two-component BECs with attractive intra-component 
interactions. General properties of vector solitons and 
their stability are studied variationally and numerically 
for both attractive and repulsive inter-component in- 
teractions. We found different types of soliton-soliton 
pairs including phase-separated pairs where one com- 
ponent is pushed outwards and forms a ring-like shell 
and the other component is compressed due to repul- 
sive inter-component interactions. It turns out that for 
some values of the chemical potentials phase-separated 
steady-states coexist with collocated states characterized 
by bell-shaped density distributions in both components. 

Wc performed a linear stability analysis of small az- 
imuthal perturbations and checked these results by an 
extensive series of numerical simulations. For attractive 
inter-component interactions matter-wave bright vector 
solitons are demonstrated to be stable throughout the ex- 
istence domain. For BEC components which repel each 
other various unstable evolution scenarios including col- 
lapse and azimuthal symmetry-breaking instabilities are 
observed. The instabilities, as a rule, lead either to sep- 
aration of the condensed phases and then a collapse of 
the stronger over-critical {N > Ncr) component or a pe- 
riodic relative motion of the under-critical (N < Ncr) 
solitonic components backwards and forwards near the 
bottom of trapping potential. Nevertheless, there are 
conditions where complete stabilization of vector solitons 
is observed even in the case of repulsive intcr-component 
interactions. 
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